Differential Binding Analysis of ChIP Seq data GSE95632
Avantika Diwadkar (avantika.diwadkar@pennmedicine.upenn.edu). Using code chunks from DiffBind.Rmd by Mengyuan Kan (mengykan@upenn.edu)
SETUP
Installing the following packages if not already present
LOADING LIBRARIES
Loading the required libraries
Read in Phenotype file with SRR numbers
pheno_file = read.table("results/GSE95632_Phenotype_withoutQC.txt", sep="\t")
colnames(pheno_file) <- c("SRR","SampleID","ID","Species","Tissue","Treatment","genome","Antibody","Ab+Treatment")
DT::datatable(pheno_file, rownames=FALSE, options = list(columnDefs = list(list(className = 'dt-center', targets = "_all"))))Subset samples and add other columns to make .csv file
#Remove input controls
pheno.main <- pheno_file[-c(2,5,17,18),] %>%select("SampleID")
#Peak file paths
Ab = as.factor(gsub(" .*","",pheno.main$SampleID))
Code = as.factor(gsub(".*_","",pheno.main$SampleID))
s_names <- as.factor(paste(Ab,Code,sep = "_"))
#Wrangle and edit columns to get necessary format
#SampleID
pheno.main$SampleID <- as.factor(gsub(".*_","",pheno.main$SampleID))
#Tissue
pheno.main$Tissue <- as.factor("ASM")
#Factor
pheno.main$Factor<- as.factor(pheno_file$Antibody[-c(2,5,17,18)])
#Treatment
pheno.main$Condition <- as.factor(c("control","control","dex","dex","dex_GFP","dex_GFP","dex_KLF","dex_KLF","control","control","dex","dex","dex_GFP","dex_GFP","dex_KLF","dex_KLF"))
pheno.main$Treatment <- as.factor(pheno.main$Condition)
#Replicates
pheno.main$Replicate <- as.factor(c(1,2))
#ControlID
pheno.main$ControlID <- as.factor(c("1D_input","1D_input","2A_input","2A_input","2A_input","2A_input","2A_input","2A_input"))
#File paths
srr_nos <- pheno_file[-c(2,5,17,18),]%>%select("SRR")
pheno.main$bamReads <- as.factor(paste("/Users/diwadkar/Desktop/ChIP-Seq/GSE95632/data/bam/",srr_nos$SRR,".bam",sep=""))
pheno.main$bamControl <- as.factor(paste("/Users/diwadkar/Desktop/ChIP-Seq/GSE95632/data/bam/",c("SRR5309352","SRR5309352","SRR5309355","SRR5309355","SRR5309355","SRR5309355","SRR5309355","SRR5309355"),".bam",sep=""))
pheno.main$Peaks <- as.factor(paste("/Users/diwadkar/Desktop/ChIP-Seq/GSE95632/data/macs2/",s_names,"_peaks.narrowPeak",sep=""))
pheno.main$PeakCaller <- as.factor("bed")
rownames(pheno.main) <- seq(length=nrow(pheno.main))
#View
#GR
pheno.gr <- pheno.main %>% filter(Factor=="GR")
DT::datatable(pheno.gr, rownames=FALSE, options = list(columnDefs = list(list(className = 'dt-center', targets = "_all"))))
write.csv(pheno.gr,"~/Desktop/ChIP-Seq/GSE95632/GSE95632info_gr.csv",row.names=FALSE)
#RNAP2
pheno.rp <- pheno.main %>% filter(Factor=="RNAP2")
DT::datatable(pheno.rp, rownames=FALSE, options = list(columnDefs = list(list(className = 'dt-center', targets = "_all"))))
write.csv(pheno.rp,"~/Desktop/ChIP-Seq/GSE95632/GSE95632info_rp.csv",row.names=FALSE)Obtaining differentially bound sites
Reading in the peaksets
Set up an experiment to analyze is with a sample sheet (.csv). The peaksets are read in using the dba function.
The metadata shows how many peaks (108812 sites) after merging overlapping ones (140789 total), and the dimensions of dba.plotPCA (default binding matrix of 16 samples by the N sites that overlap in at least two of the samples).
A correlation heatmap gives an initial clustering of the samples using the cross-correlations of each row of the binding matrix.
#GR
dat.gr <- dba(sampleSheet="/Users/diwadkar/Desktop/ChIP-Seq/GSE95632/GSE95632info_gr.csv") ## 1D ASM GR control control 1 bed## 1E ASM GR control control 2 bed## 2A ASM GR dex dex 1 bed## 2E ASM GR dex dex 2 bed## 3D ASM GR dex_GFP dex_GFP 1 bed## 3E ASM GR dex_GFP dex_GFP 2 bed## 4B ASM GR dex_KLF dex_KLF 1 bed## 4C ASM GR dex_KLF dex_KLF 2 bedplot(dat.gr)
#RNAP2
dat.rp <- dba(sampleSheet="/Users/diwadkar/Desktop/ChIP-Seq/GSE95632/GSE95632info_rp.csv") ## 5E ASM RNAP2 control control 1 bed## 5b ASM RNAP2 control control 2 bed## 6d ASM RNAP2 dex dex 1 bed## 6e ASM RNAP2 dex dex 2 bed## 7c ASM RNAP2 dex_GFP dex_GFP 1 bed## 7e ASM RNAP2 dex_GFP dex_GFP 2 bed## 9b ASM RNAP2 dex_KLF dex_KLF 1 bed## 9e ASM RNAP2 dex_KLF dex_KLF 2 bedplot(dat.rp)
Counting reads
The next step is to calculate a binding matrix with scores based on read counts for every sample (affinity scores), rather than confidence scores for only those peaks called in a specific sample (occupancy scores).
Use dba.count function. The current study is based on a transcription factor GR that binds to the DNA, resulting in “punctate”, narrow peaks, it is advisable to use the “summits” option to re-center each peak around the point of greatest enrichment. This keeps the peaks at a consistent width (in this case, with summits=250, the peaks will be 500bp, extending 250bp up- and down- stream of the summit).
If you do not have the raw reads available to you, this object is available loading using e.g. data(tamoxifen_counts).
#GR
dat.gr <- dba.count(dat.gr,summits=250) # This step takes a while## Re-centering peaks...#RNAP2
dat.rp <- dba.count(dat.rp,summits=250)## Re-centering peaks...This shows that all the samples are using the same, 82032 length consensus peakset. Also, a new column has been added, called FRiP (Fraction of Reads in Peaks). This is the proportion of reads for that sample that overlap a peak in the consensus peakset, and can be used to indicate which samples show more enrichment overall.
dat.gr## 8 Samples, 98680 sites in matrix:
## ID Tissue Factor Condition Treatment Replicate Caller Intervals FRiP
## 1 1D ASM GR control control 1 counts 98680 0.30
## 2 1E ASM GR control control 2 counts 98680 0.30
## 3 2A ASM GR dex dex 1 counts 98680 0.31
## 4 2E ASM GR dex dex 2 counts 98680 0.35
## 5 3D ASM GR dex_GFP dex_GFP 1 counts 98680 0.29
## 6 3E ASM GR dex_GFP dex_GFP 2 counts 98680 0.28
## 7 4B ASM GR dex_KLF dex_KLF 1 counts 98680 0.28
## 8 4C ASM GR dex_KLF dex_KLF 2 counts 98680 0.25dat.rp## 8 Samples, 29852 sites in matrix:
## ID Tissue Factor Condition Treatment Replicate Caller Intervals FRiP
## 1 5E ASM RNAP2 control control 1 counts 29852 0.09
## 2 5b ASM RNAP2 control control 2 counts 29852 0.10
## 3 6d ASM RNAP2 dex dex 1 counts 29852 0.11
## 4 6e ASM RNAP2 dex dex 2 counts 29852 0.10
## 5 7c ASM RNAP2 dex_GFP dex_GFP 1 counts 29852 0.09
## 6 7e ASM RNAP2 dex_GFP dex_GFP 2 counts 29852 0.06
## 7 9b ASM RNAP2 dex_KLF dex_KLF 1 counts 29852 0.08
## 8 9e ASM RNAP2 dex_KLF dex_KLF 2 counts 29852 0.09Plot a new correlation heatmap based on the affinity scores.
plot(dat.gr)
plot(dat.rp)
Establishing a contrast
Before running the differential analysis, we need to tell DiffBind which cell lines fall in which groups.
This uses the Condition metadata (dex vs. control) to set up a contrast with 2 samples in the dex group and 2 samples in control group.
The comparison indicates the group 1 (dex) vs. group 2 (control)
#GR
dat.gr <- dba.contrast(dat.gr, group1=dat.gr$masks$dex, group2=dat.gr$masks$control, name1="dex", name2="control")
dat.gr <- dba.contrast(dat.gr, group1=dat.gr$masks$dex_GFP, group2=dat.gr$masks$dex_KLF, name1="dex_GFP", name2="dex_KLF")
#dat <- dba.contrast(dat, categories=DBA_CONDITION, minMembers=2)
#RNAP2
dat.rp <- dba.contrast(dat.rp, group1=dat.rp$masks$dex, group2=dat.rp$masks$control, name1="dex", name2="control")
dat.rp <- dba.contrast(dat.rp, group1=dat.rp$masks$dex_GFP, group2=dat.rp$masks$dex_KLF, name1="dex_GFP", name2="dex_KLF")
#dat <- dba.contrast(dat, categories=DBA_CONDITION, minMembers=2)Performing the differential analysis
The default analysis will run an DESeq2 analysis using the default binding matrix. The resultant DBA object shows that 37582 of the 82032 sites are identified as being significantly differentially bound (DB) using a FDR <= 0.05. A correlation heatmap can be plotted, based on the result of the analysis.
#GR - dex vs control
dat.gr <- dba.analyze(dat.gr)
dat.gr## 8 Samples, 98680 sites in matrix:
## ID Tissue Factor Condition Treatment Replicate Caller Intervals FRiP
## 1 1D ASM GR control control 1 counts 98680 0.30
## 2 1E ASM GR control control 2 counts 98680 0.30
## 3 2A ASM GR dex dex 1 counts 98680 0.31
## 4 2E ASM GR dex dex 2 counts 98680 0.35
## 5 3D ASM GR dex_GFP dex_GFP 1 counts 98680 0.29
## 6 3E ASM GR dex_GFP dex_GFP 2 counts 98680 0.28
## 7 4B ASM GR dex_KLF dex_KLF 1 counts 98680 0.28
## 8 4C ASM GR dex_KLF dex_KLF 2 counts 98680 0.25
##
## 2 Contrasts:
## Group1 Members1 Group2 Members2 DB.DESeq2
## 1 dex 2 control 2 39099
## 2 dex_GFP 2 dex_KLF 2 28861plot(dat.gr, contrast=1)
plot(dat.gr, contrast=2)
#GR - GFP dex vs KLF dex
#dat.gr2 <- dba.analyze(dat.gr2)
#dat.gr2
#plot(dat.gr2, contrast=1)
#RNAP2 - dex vs control
dat.rp <- dba.analyze(dat.rp)
dat.rp## 8 Samples, 29852 sites in matrix:
## ID Tissue Factor Condition Treatment Replicate Caller Intervals FRiP
## 1 5E ASM RNAP2 control control 1 counts 29852 0.09
## 2 5b ASM RNAP2 control control 2 counts 29852 0.10
## 3 6d ASM RNAP2 dex dex 1 counts 29852 0.11
## 4 6e ASM RNAP2 dex dex 2 counts 29852 0.10
## 5 7c ASM RNAP2 dex_GFP dex_GFP 1 counts 29852 0.09
## 6 7e ASM RNAP2 dex_GFP dex_GFP 2 counts 29852 0.06
## 7 9b ASM RNAP2 dex_KLF dex_KLF 1 counts 29852 0.08
## 8 9e ASM RNAP2 dex_KLF dex_KLF 2 counts 29852 0.09
##
## 2 Contrasts:
## Group1 Members1 Group2 Members2 DB.DESeq2
## 1 dex 2 control 2 1616
## 2 dex_GFP 2 dex_KLF 2 235plot(dat.rp, contrast=1)
plot(dat.rp, contrast=2)
#RNAP2 - GFP dex vs KLF dex
#dat.rp2 <- dba.analyze(dat.rp2)
#dat.rp2
#plot(dat.rp2, contrast=1)Retrieve the differentially bound sites
The output is a GRanges object, appropriate for downstream processing.
The value columns show:
- the mean read concentration over all the samples: the default calculation uses log2 normalized ChIP read counts with control read counts subtracted
- the mean concentration of the first (dex) group
- the mean concentration of second (control) group
- the Fold column shows the difference in mean concentrations between the two groups (Conc_dex - Conc_control): a positive value indicating increased binding affinity in the dex group and a negative value indicating increased binding affinity in the control group.
- a raw p-value: confidence measures for identifying these sites as differentially bound
- FDR: a multiple testing corrected FDR in the final column.
#GR1
dat.gr.DB <- dba.report(dat.gr)
dat.gr.DB## GRanges object with 39099 ranges and 6 metadata columns:
## seqnames ranges strand | Conc Conc_dex
## <Rle> <IRanges> <Rle> | <numeric> <numeric>
## 17705 chr11 114179230-114179730 * | 9.44 10.41
## 11837 chr10 74056322-74056822 * | 9.12 10.06
## 42623 chr2 20438541-20439041 * | 8.75 9.67
## 17699 chr11 114154432-114154932 * | 8.78 9.72
## 36799 chr17 59825482-59825982 * | 8.69 9.61
## ... ... ... ... . ... ...
## 8163 chr1 221989640-221990140 * | 5.86 5.4
## 73147 chr5 137188108-137188608 * | 2.83 3.61
## 88784 chr8 67245151-67245651 * | 2.83 3.61
## 11693 chr10 71230593-71231093 * | 6.1 6.41
## 25151 chr13 74229142-74229642 * | 3.35 4.04
## Conc_control Fold p-value FDR
## <numeric> <numeric> <numeric> <numeric>
## 17705 4.71 5.7 1.16e-119 1.15e-114
## 11837 5.36 4.7 2.93e-111 1.44e-106
## 42623 5.38 4.29 6.53e-93 2.15e-88
## 17699 5.26 4.46 4.36e-92 1.08e-87
## 36799 5.51 4.1 7.36e-83 1.45e-78
## ... ... ... ... ...
## 8163 6.2 -0.8 0.0198 0.05
## 73147 0.98 2.63 0.0198 0.05
## 88784 0.98 2.63 0.0198 0.05
## 11693 5.69 0.72 0.0198 0.05
## 25151 1.96 2.08 0.0198 0.05
## -------
## seqinfo: 33 sequences from an unspecified genome; no seqlengths#GR2
#dat.gr2.DB <- dba.report(dat.gr2)
#dat.gr2.DB
#RNAP2
dat.rp.DB <- dba.report(dat.rp)
dat.rp.DB## GRanges object with 1616 ranges and 6 metadata columns:
## seqnames ranges strand | Conc Conc_dex
## <Rle> <IRanges> <Rle> | <numeric> <numeric>
## 17681 chr20 47318285-47318785 * | 7.84 8.63
## 10608 chr16 56608395-56608895 * | 9.2 9.79
## 18463 chr22 30246411-30246911 * | 7.36 5.79
## 23916 chr6 35727765-35728265 * | 6.22 7.1
## 28928 chr9 129320831-129321331 * | 7.59 8.26
## ... ... ... ... . ... ...
## 23416 chr5 179829909-179830409 * | 4.26 3.07
## 10695 chr16 66996370-66996870 * | 3.43 4.23
## 3982 chr10 110535930-110536430 * | 3.41 1.32
## 17476 chr20 35443015-35443515 * | 3.59 4.4
## 20292 chr3 155079695-155080195 * | 4.16 4.83
## Conc_control Fold p-value FDR
## <numeric> <numeric> <numeric> <numeric>
## 17681 5.98 2.64 1.95e-36 5.47e-32
## 10608 8.17 1.63 2.52e-25 3.55e-21
## 18463 8.1 -2.3 4.13e-24 3.87e-20
## 23916 3.65 3.45 1.94e-19 1.36e-15
## 28928 6.29 1.97 2.59e-18 1.46e-14
## ... ... ... ... ...
## 23416 4.9 -1.83 0.00286 0.0499
## 10695 1.43 2.8 0.00286 0.0499
## 3982 4.23 -2.91 0.00287 0.0499
## 17476 1.53 2.87 0.00287 0.0499
## 20292 2.84 1.99 0.00287 0.0499
## -------
## seqinfo: 24 sequences from an unspecified genome; no seqlengths#RNAP2
#dat.rp2.DB <- dba.report(dat.rp2)
#dat.rp2.DBVisualization
PCA
A PCA plot shows a better understanding of clustering and association between samples than heatmaps.
#PCA plot for normalized read counts for all binding sites
dba.plotPCA(dat.gr,DBA_TREATMENT,label=DBA_CONDITION)
#dba.plotPCA(dat.gr2,DBA_TREATMENT,label=DBA_CONDITION)
#dba.plotPCA(dat,DBA_FACTOR,label=DBA_CONDITION)
#dba.plotPCA(dat.rp1,DBA_TREATMENT,label=DBA_CONDITION)
dba.plotPCA(dat.rp,DBA_TREATMENT,label=DBA_CONDITION)
#PCA plot using only the differentially bound sites, using an FDR threshold of 0.05
dba.plotPCA(dat.gr,contrast = 1,label=DBA_REPLICATE)
dba.plotPCA(dat.gr,contrast = 2,label=DBA_REPLICATE)
#dba.plotPCA(dat.gr1,contrast = 1,label=DBA_TREATMENT)
#dba.plotPCA(dat.gr2,contrast = 1,label=DBA_REPLICATE)
#dba.plotPCA(dat.gr2,contrast = 1,label=DBA_TREATMENT)
#dba.plotPCA(dat,contrast = 1,label=DBA_FACTOR)
#RNAP2
dba.plotPCA(dat.rp,contrast = 1,label=DBA_REPLICATE)
dba.plotPCA(dat.rp,contrast = 2,label=DBA_REPLICATE)
#dba.plotPCA(dat.rp1,contrast = 1,label=DBA_TREATMENT)
#dba.plotPCA(dat.rp2,contrast = 1,label=DBA_REPLICATE)
#dba.plotPCA(dat.rp2,contrast = 1,label=DBA_TREATMENT)MA PLOTS
MA plots are a useful way to visualize the effect of normalization on data, as well as seeing which of the datapoints are being identified as differentially bound. An MA plot can be obtained for the control vs dex treatment cell:
dba.plotMA(dat.gr)
dba.plotMA(dat.rp)
#Same data shown with the concentrations of each sample groups plotted against each other -
dba.plotMA(dat.gr, bXY=TRUE)
#dba.plotMA(dat.gr2, bXY=TRUE)
dba.plotMA(dat.rp, bXY=TRUE)
#dba.plotMA(dat.rp2, bXY=TRUE)VOLCANO PLOTS
Volcano plots also highlight significantly differentially bound sites and show their fold changes.
dba.plotVolcano(dat.gr)
#dba.plotVolcano(dat.gr2)
dba.plotVolcano(dat.rp)
#dba.plotVolcano(dat.rp2)BOXPLOTS
Box plots of read distributions for significantly differentially bound (DB) sites - The plot shows in the first two boxes that amongst differentially bound sites overall, the dex samples have a somewhat higher mean read concentration. The next two boxes show the distribution of reads in differentially bound sites that exhibit increased affinity in the dex samples, while the final two boxes show the distribution of reads in differentially bound sites that exhibit increased affinity in the Control samples.
Pvals is a matrix of p-values (computed using a two-sided Wilcoxon ‘MannWhitney’ test, paired where appropriate) indicating which of these distributions are significantly different from another distribution.
#GR1
sum(dat.gr.DB$Fold<0)## [1] 15983sum(dat.gr.DB$Fold>0)## [1] 23116pvals1 <- dba.plotBox(dat.gr)
pvals1## dex.DB control.DB dex.DB+ control.DB+ dex.DB-
## dex.DB 1.00e+00 0.00e+00 2.82e-212 0 9.73e-131
## control.DB 0.00e+00 1.00e+00 4.24e-120 0 0.00e+00
## dex.DB+ 2.82e-212 4.24e-120 1.00e+00 0 0.00e+00
## control.DB+ 0.00e+00 0.00e+00 0.00e+00 1 0.00e+00
## dex.DB- 9.73e-131 0.00e+00 0.00e+00 0 1.00e+00
## control.DB- 0.00e+00 0.00e+00 0.00e+00 0 0.00e+00
## control.DB-
## dex.DB 0
## control.DB 0
## dex.DB+ 0
## control.DB+ 0
## dex.DB- 0
## control.DB- 1#GR2
#sum(dat.gr2.DB$Fold<0)
#sum(dat.gr2.DB$Fold>0)
#pvals2 <- dba.plotBox(dat.gr2)
#pvals2
#RNAP2-1
sum(dat.rp.DB$Fold<0)## [1] 576sum(dat.rp.DB$Fold>0)## [1] 1040pvals3 <- dba.plotBox(dat.rp)
pvals3## dex.DB control.DB dex.DB+ control.DB+ dex.DB-
## dex.DB 1.00e+00 2.58e-33 1.21e-74 7.23e-22 4.39e-35
## control.DB 2.58e-33 1.00e+00 4.67e-17 1.86e-52 2.65e-83
## dex.DB+ 1.21e-74 4.67e-17 1.00e+00 4.96e-96 3.85e-155
## control.DB+ 7.23e-22 1.86e-52 4.96e-96 1.00e+00 4.97e-01
## dex.DB- 4.39e-35 2.65e-83 3.85e-155 4.97e-01 1.00e+00
## control.DB- 2.12e-90 5.44e-25 1.41e-02 1.18e-108 9.99e-172
## control.DB-
## dex.DB 2.12e-90
## control.DB 5.44e-25
## dex.DB+ 1.41e-02
## control.DB+ 1.18e-108
## dex.DB- 9.99e-172
## control.DB- 1.00e+00#RNAP2-2
#sum(dat.rp2.DB$Fold<0)
#sum(dat.rp2.DB$Fold>0)
#pvals4 <- dba.plotBox(dat.rp2)
#pvals4HEATMAPS
#GR
#Simple heatmap
corvals1 <- dba.plotHeatmap(dat.gr)
#Binding affinity heatmap
corvals1 <- dba.plotHeatmap(dat.gr, contrast=1, correlations=FALSE)
corvals1 <- dba.plotHeatmap(dat.gr, contrast=2, correlations=FALSE)
#Simple heatmap
#corvals2 <- dba.plotHeatmap(dat.gr2)
#Binding affinity heatmap
#corvals2 <- dba.plotHeatmap(dat.gr2, contrast=1, correlations=FALSE)
#RNAP2
#Simple heatmap
corvals3 <- dba.plotHeatmap(dat.rp)
#Binding affinity heatmap
corvals3 <- dba.plotHeatmap(dat.rp, contrast=1, correlations=FALSE)
corvals3 <- dba.plotHeatmap(dat.rp, contrast=2, correlations=FALSE)
#Simple heatmap
#corvals4 <- dba.plotHeatmap(dat.rp2)
#Binding affinity heatmap
#corvals4 <- dba.plotHeatmap(dat.rp2, contrast=1, correlations=FALSE)Session information
R version 3.5.0 (2018-04-23)
**Platform:** x86_64-apple-darwin15.6.0 (64-bit)
locale: en_US.UTF-8||en_US.UTF-8||en_US.UTF-8||C||en_US.UTF-8||en_US.UTF-8
attached base packages: parallel, stats4, stats, graphics, grDevices, utils, datasets, methods and base
other attached packages: bindrcpp(v.0.2.2), DT(v.0.4), dplyr(v.0.7.6), plyr(v.1.8.4), pander(v.0.6.1), DiffBind(v.2.8.0), SummarizedExperiment(v.1.10.1), DelayedArray(v.0.6.1), BiocParallel(v.1.14.1), matrixStats(v.0.53.1), Biobase(v.2.40.0), GenomicRanges(v.1.32.3), GenomeInfoDb(v.1.16.0), IRanges(v.2.14.10), S4Vectors(v.0.18.3) and BiocGenerics(v.0.26.0)
loaded via a namespace (and not attached): amap(v.0.8-16), colorspace(v.1.3-2), rjson(v.0.2.20), hwriter(v.1.3.2), rprojroot(v.1.3-2), htmlTable(v.1.12), XVector(v.0.20.0), base64enc(v.0.1-3), rstudioapi(v.0.7), ggrepel(v.0.8.0), bit64(v.0.9-7), AnnotationDbi(v.1.42.1), splines(v.3.5.0), geneplotter(v.1.58.0), knitr(v.1.20), Formula(v.1.2-3), Rsamtools(v.1.32.0), annotate(v.1.58.0), cluster(v.2.0.7-1), GO.db(v.3.6.0), pheatmap(v.1.0.10), graph(v.1.58.0), compiler(v.3.5.0), httr(v.1.3.1), GOstats(v.2.46.0), backports(v.1.1.2), assertthat(v.0.2.0), Matrix(v.1.2-14), lazyeval(v.0.2.1), limma(v.3.36.2), acepack(v.1.4.1), htmltools(v.0.3.6), prettyunits(v.1.0.2), tools(v.3.5.0), gtable(v.0.2.0), glue(v.1.2.0), GenomeInfoDbData(v.1.1.0), Category(v.2.46.0), systemPipeR(v.1.14.0), ShortRead(v.1.38.0), Rcpp(v.0.12.17), Biostrings(v.2.48.0), gdata(v.2.18.0), rtracklayer(v.1.40.3), stringr(v.1.3.1), gtools(v.3.8.1), XML(v.3.98-1.11), edgeR(v.3.22.3), zlibbioc(v.1.26.0), scales(v.0.5.0), hms(v.0.4.2), RBGL(v.1.56.0), RColorBrewer(v.1.1-2), BBmisc(v.1.11), yaml(v.2.1.19), gridExtra(v.2.3), memoise(v.1.1.0), ggplot2(v.2.2.1), rpart(v.4.1-13), biomaRt(v.2.36.1), latticeExtra(v.0.6-28), stringi(v.1.2.3), RSQLite(v.2.1.1), genefilter(v.1.62.0), checkmate(v.1.8.5), GenomicFeatures(v.1.32.0), caTools(v.1.17.1), rlang(v.0.2.1), pkgconfig(v.2.0.1), BatchJobs(v.1.7), bitops(v.1.0-6), evaluate(v.0.10.1), lattice(v.0.20-35), purrr(v.0.2.5), bindr(v.0.1.1), labeling(v.0.3), GenomicAlignments(v.1.16.0), htmlwidgets(v.1.2), bit(v.1.1-14), tidyselect(v.0.2.4), GSEABase(v.1.42.0), AnnotationForge(v.1.22.0), magrittr(v.1.5), sendmailR(v.1.2-1), DESeq2(v.1.20.0), R6(v.2.2.2), gplots(v.3.0.1), Hmisc(v.4.1-1), DBI(v.1.0.0), foreign(v.0.8-70), pillar(v.1.2.3), nnet(v.7.3-12), survival(v.2.42-4), RCurl(v.1.95-4.10), tibble(v.1.4.2), crayon(v.1.3.4), KernSmooth(v.2.23-15), rmarkdown(v.1.10), progress(v.1.2.0), locfit(v.1.5-9.1), grid(v.3.5.0), data.table(v.1.11.4), blob(v.1.1.1), Rgraphviz(v.2.24.0), digest(v.0.6.15), xtable(v.1.8-2), brew(v.1.0-6) and munsell(v.0.5.0)